行列の最小二乗法について調べてみた
カフェチームの山本です。
現在、カフェチームでは、機械学習を用いて、画像内に映っている人物の骨格や手を検出することで、どの商品を取り出したかを判定する方法を検討しています。
前回は、実際にカメラの位置や回転を計測することなく、この外部パラメータを推定する方法として、最小二乗法を用いました。
今回は、上記の記事を作成するにあたって調べた最小二乗法について、計算過程や証明方法をまとめました。結論としては、どのような次元でも、正規方程式を用いて計算できることがわかりました。
お詫び
本ページの数式が正しく表示されていません。現在、修正中です。数式まで見たい方はこちらのページをご覧ください。
０．最小二乗法の概要
説明
例として、２次元平面において、[(x0, y0), (x1, y1), ...]のような対応関係がある、サンプル（測定点）があるとします。このサンプルを直線で近似することを考えます。すなわち、以下のようなモデルです。
[latex]y=qx+p[/latex]
この近似において、「各測定点でのy座標の誤差の２乗」の和を最も小さくするような(q, p)を求めるのが最小二乗法です。これによって、測定点の近くを通る直線を求めることができます。
１．最小二乗法の定式化・計算
（１）入力：１元・出力：１元の場合
ステップ０：誤差の定式化
誤差２乗の和を定式化すると、以下のように表すことができます。
[latex]E = \left\{y_0-(qx_0+p)\right\}^2+\left\{y_1-(qx_1+p)\right\}^2+\cdots[/latex]
ステップ１：行列で書き直す
これは、行列を使って以下のように書き直すことができます。
[latex]b = \left( \begin{matrix} y_0 \\ y_1 \\ \vdots \end{matrix} \right)[/latex]
[latex]A = \left( \begin{matrix} x_0 & 1 \\ x_1 & 1 \\ \vdots \end{matrix} \right)[/latex]
[latex]x = \left( \begin{matrix} q \\ p \end{matrix} \right)[/latex]
[latex]E = \left\| b - Ax \right\|^2[/latex]
（このEは、下に示すように上の式と同じです）
[latex]\begin{matrix} E &=& \left\| b - Ax \right\|^2 \\ &=& \left\| \left(\begin{matrix} y_0 - (qx_0+p) \\ y_1 - (qx_1+p) \\ \vdots \end{matrix}\right)\right\|^2 \\ &=& \left\{y_0-(qx_0+p)\right\}^2+\left\{y_1-(qx_1+p)\right\}^2+\cdots \end{matrix}[/latex]
ステップ２：Eを変形させる
Eを変形させると以下のようになります。
[latex]\begin{matrix} E &=& \left\| b - Ax \right\|^2 \\ &=& \left(b - Ax \right)^T \left(b - Ax \right) \\ &=& \left(b^T - (Ax)^T \right) \left(b - Ax \right) \\ &=& \left(b^T - x^TA^T \right)\left(b - Ax \right) \\ &=& b^Tb- x^TA^Tb - b^TAx + x^TA^TAx \end{matrix}[/latex]
- 補足１：n行1列ベクトルx, yがあるとき、以下の関係があります。
[latex]x^Ty = \sum_ix_iy_i=y^Tx[/latex]
そのため、上のEの変形における最後の式の2項目と3項目は同じ値です。よって下のように変形できます。
[latex]\begin{matrix} E &=& b^Tb- x^TA^Tb - b^TAx + x^TA^TAx \\ &=& b^Tb- 2b^TAx + x^TA^TAx \end{matrix}[/latex]
これを各成分を代入して表すと以下のようになります。
[latex]\begin{matrix} E &=& b^Tb- x^TA^Tb - b^TAx + x^TA^TAx \\ &=& \sum_i b_i^2 - \sum_i y_i(qx_i+p) - \sum_i y_i(qx_i+p) + \sum_i(qx_i+p)^2 \\ &=& \sum_i b_i^2 - 2\sum_i y_i(qx_i+p) + \sum_i(qx_i+p)^2 \\ &=& \sum_i b_i^2 - 2\sum_i y_ix_iq - 2\sum_i y_ip + \sum_ix_i^2q^2 + 2\sum_i x_i qp + \sum_ip^2 \end{matrix}[/latex]
ステップ３：偏微分・極小点（成分計算の場合）
上の式から、Eはpの2次式であり、かつ、qの2次式であることがわかります。また、pとqの2次の項が係数が正です。そのため、Eをq,pそれぞれで偏微分した値が同時に0になる点があれば、その点でEが最小になることがわかります。
[latex]\left(\begin{matrix} \frac{\partial E}{\partial q} \\ \frac{\partial E}{\partial p} \end{matrix}\right)=0[/latex]
[latex]\Leftrightarrow \left(\begin{matrix} - 2\sum_iy_ix_i + \sum_i x_i^2 2q + 2 \sum_i x_i p \\ - 2 \sum_i y_i + 2\sum_i x_i q + \sum_i 2p \end{matrix} \right)=0[/latex]
[latex]\Leftrightarrow \left(\begin{matrix} - 2 \sum_iy_ix_i + 2 \sum_i x_i^2 q + 2 \sum_i x_i p \\ - 2 \sum_i y_i + 2 \sum_i x_i q + 2 \sum_i p \end{matrix} \right)=0[/latex]
[latex]\Leftrightarrow \left(\begin{matrix} - \sum_iy_ix_i + \sum_i x_i^2 q + \sum_i x_i p \\ - \sum_i y_i + \sum_i x_i q + \sum_i p \end{matrix} \right)=0[/latex]
[latex]\Leftrightarrow \left( \begin{matrix} \sum_i x_i^2 q + \sum_i x_i p \\ \sum_i x_i q + \sum_i p \end{matrix} \right) = \left( \begin{matrix} \sum_iy_ix_i \\ \sum_i y_i \end{matrix} \right)[/latex]
[latex]\Leftrightarrow \left( \begin{matrix} \sum_i x_i^2 & \sum_i x_i \\ \sum_i x_i & \sum_i 1 \end{matrix} \right) \left( \begin{matrix} q \\ p\end{matrix} \right) = \left( \begin{matrix} \sum_iy_ix_i \\ \sum_i y_i \end{matrix} \right)[/latex]
左辺の第一項が正則行列である場合、
[latex]\Leftrightarrow \left( \begin{matrix} q \\ p\end{matrix} \right) = \left( \begin{matrix} \sum_i x_i^2 & \sum_i x_i \\ \sum_i x_i & \sum_i 1 \end{matrix} \right)^{-1}\left( \begin{matrix} \sum_iy_ix_i \\ \sum_i y_i \end{matrix} \right)[/latex]
上のように式変形することで、(q,p)を求めることができます。
ステップ３：偏微分・極小点（行列計算の場合）
これをxの各成分に関して偏微分すると以下のようになります。
[latex]\nabla = \left( \begin{matrix} \frac{\partial}{\partial q} \\ \frac{\partial}{\partial p} \end{matrix} \right)[/latex]
[latex]\begin{matrix} \nabla E &=& - 2 \left( b^T A\right)^T + \left( A^T A + (A^TA)^T \right)x \\ &=& - 2 A^T b + \left( A^T A + A^T A \right)x \\ &=& - 2 A^T b + 2 A^T A x \end{matrix}[/latex]
Eが最小になるには、偏微分の値が0になることが必要条件です。
[latex]\nabla E = 0[/latex]
[latex]\Leftrightarrow - 2 A^T b + 2 A^T A x = 0[/latex]
[latex]\Leftrightarrow 2 A^T A x = 2 A^T b[/latex]
[latex]\Leftrightarrow A^T A x = A^T b[/latex]
という式が導けます。A^{T} Aが正則行列のとき、
[latex]x = \left( A^T A \right)^{-1} A^T b[/latex]
として、xを計算することができます。これは、上の成分を代入して計算した場合と同じです。
- 補足３：A^{T}Aが正則行列がどうかの判定
A^{T}Aが正則行列であることは、以下のように考えられます。
正則行列の性質から、以下が成り立ちます。
[latex]A^{T}Aが正則である\Leftrightarrow \rm{rank}(A^{T}A)=2[/latex]
ランク（階数）の性質から、以下が成り立ちます
[latex]\rm{rank}(A^{T}A)=2 \Leftrightarrow \rm{rank}(A)=2[/latex]
また、階数の定義から、以下が成り立ちます。
[latex]\rm{rank}(A)=2 \Leftrightarrow Aの行ベクトルにおいて、線形独立なものが２つ以上存在する[/latex]
Aは計測したサンプルのx座標を並べた行列です。なので、上と組み合わせると、x座標の異なる点で２つ以上サンプルを計測すれば、A^{T}Aが正則行列になることがわかります。
（つまり、２点以上測定して直線を引けるようにすれば良いということです）
（２）入力：多元・出力：１元の場合
前節は入力が１元（ｘ）、出力が１元（ｙ）でしたが、入力が多元になっても同様に計算できます。
以下では、入力をx,yとして、出力をzとします
誤差２乗の和を定式化すると、以下のようです。
[latex]E = \left\{z_0-(ry_0 + qx_0+p)\right\}^2+\left\{z_1-(ry_1+qx_1+p)\right\}^2+\cdots[/latex]
これは、行列を使って以下のように書き直すことができます。
[latex]b = \left( \begin{matrix} z_0 \\ z_1 \\ \vdots \end{matrix} \right)[/latex]
[latex]A = \left( \begin{matrix} y_0 & x_0 & 1 \\ y_1 & x_1 & 1 \\ \vdots \end{matrix} \right)[/latex]
[latex]x = \left( \begin{matrix} r \\ q \\ p \end{matrix} \right)[/latex]
とすると、Eは以下のようになり、前節と同じ形になります。
[latex]E = \left\| b - Ax \right\|^2[/latex]
よって、Eを最小にするxは、同様に計算できます。
[latex]x = \left( A^T A \right)^{-1} A^T b[/latex]
（３）入力：多元・出力：多元の場合
出力が多元になっても同様の計算ができます。
以下では、入力を(x,y,z)、出力を(x',y',z')とします。
誤差２乗の和を定式化すると、以下のようです。
[latex]E = \left\{z_0'-(r_{11}z_0 + r_{21}y_0 + r_{31}x_0 +r_{41})\right\}^2 + \left\{y_0'-(r_{12}z_0 + r_{22}y_0 + r_{32}x_0 + r_{42})\right\}^2 + \left\{x_0'-(r_{13}z_0 + r_{23}y_0 + r_{33}x_0 + r_{43})\right\}^2 + \left\{z_1'-(r_{11}z_1 + r_{21}y_1 + r_{31}x_1 +r_{41})\right\}^2 + \left\{y_1'-(r_{12}z_1 + r_{22}y_1 + r_{32}x_1 + r_{42})\right\}^2 + \left\{x_1'-(r_{13}z_1 + r_{23}y_1 + r_{33}x_1 + r_{43})\right\}^2 + \cdots[/latex]
これは、行列を使って以下のように書き直すことができます。
[latex]B = \left( \begin{matrix} x_0' & y_0' & z_0' \\ x_1' & y_1' & z_1' \\ & \vdots & \end{matrix} \right) = \left( \begin{matrix} b_x & b_y & b_z \end{matrix} \right)[/latex]
[latex]A = \left( \begin{matrix} z_0 & y_0 & x_0 & 1 \\ z_1 & y_1 & x_1 & 1 \\ & \vdots & & \end{matrix} \right)[/latex]
[latex]X = \left( \begin{matrix} r_{11} & r_{12} & r_{13} \\ r_{21} & r_{22} & r_{23} \\ r_{31} & r_{32} & r_{33} \\ r_{41} & r_{42} & r_{43} \end{matrix} \right)= \left( \begin{matrix} r_1 & r_2 & r_3 \end{matrix} \right)[/latex]
このとき、Eをx,y,zそれぞれに関する誤差でまとめると、以下のようになります。
[latex]E = \left\| b_x - A r_1 \right\|^2 + \left\| b_y - A r_2 \right\|^2 + \left\| b_z - A r_3 \right\|^2[/latex]
Eの形をみると、それぞれの項でr1,r2,r3が独立しており、かつ、各項を最小化されればEが最小化されます。よって、（今回の目的である）Eを最小化するr1,r2,r3を求めるには、各項を最小化するr1,r2,r3をそれぞれ求めれば良いことになります。そしてこれは、前節と同じ計算です。
[latex]r_1 = \left( A^T A \right)^{-1} A^T b_x[/latex]
[latex]r_2 = \left( A^T A \right)^{-1} A^T b_y[/latex]
[latex]r_3 = \left( A^T A \right)^{-1} A^T b_z[/latex]
３つの式を横にならべると
[latex]\left( \begin{matrix} r_1 & r_2 & r_3 \end{matrix} \right) = \left( A^T A \right)^{-1} A^T \left( \begin{matrix} b_x & b_y & b_z \end{matrix} \right)[/latex]
[latex]\Leftrightarrow X= \left(A^T A \right)^{-1} A^T B[/latex]
となり、前節と同様の形になります。
まとめ
サンプルとなる対応関係が得られたとき、(線形な)モデルでそれらを近似する方法である、最小二乗法についてまとめました。最小二乗法では、入力出力が多次元の場合でも利用することができます。具体的には、対応関係の行列をA,Bとしたとき、以下のように計算できます。
[latex]x = \left( A^T A \right)^{-1} A^T b[/latex]
[latex]X= \left(A^T A \right)^{-1} A^T B[/latex]
参考にさせていただいたページ
最小二乗法の基礎をまとめる - エンジニアを目指す浪人のブログ
行列のランクの定義・例・性質/公式 (証明付) - 理数アラカルト -
補足
これを利用して、カメラの外部パラメータの推定を行った様子を、以下の記事で書いています。
